Curso 7 Regression Models
Curso 7 Regression Models
- Introducción
- Concepto de Regresión lineal y mínimos cuadrados
- Concepto de Residuos
- Concepto de estimación por mínimos cuadrados
- TEST 1
- Aplicación de regresión lineal y mínimos cuadrados
- Varianza de Residuos
- Inferencia
- TEST2
- Regresión multivariable
- Ajustes y relacción entre variables
- Residuos y Diagnósticos de outliers
- Selección del Modelo (VIF y ANOVA)
- Overfitting and Underfitting (ANOVA)
- TEST3
- GLM (Generalized Linear Models)
- Resultados binarios (Bernuilli)
- Resultados de conteo (Poisson)
- Ajustes en funciones
- TEST4
Introducción
En estadística, el análisis de la regresión es un proceso estadístico para estimar las relaciones entre variables. Incluye muchas técnicas para el modelado y análisis de diversas variables, cuando la atención se centra en la relación entre una variable dependiente (o Outcome o Variable de respuesta) y una o más variables independientes (o Covariables o Variables regresoras ).
Este curso cubre análisis de regresión, mínimos cuadrados e inferencia usando modelos de regresión. Los casos especiales de modelos de regresión ANOVA y ANCOVA también serán cubiertos. El análisis de los residuos y la variablidad también se investigará. Los modelos de regresión producen modelos muy facilmente interpretables, al contrario que los algoritmos de aprendizaje automático que muchas veces sacrifican interpretabilidad por un mejor rendimiento en la predicción. Esto por supuesto es el fin de toda estimación, sin embargo por simpleza e interpretabilidad deberíamos optar siempre en primer lugar por estos modelos de regresión.
Podemos clasificar los tipos de regresión según diversos criterios.
En primer lugar, en función del número de variables independientes:
Regresión simple: Cuando la variable Y depende únicamente de una única variable X.
Regresión múltiple: Cuando la variable Y depende de varias variables (X1, X2, ..., Xr)
En segundo lugar, en función del tipo de función f(X):
Regresión lineal: Cuando f(X) es una función lineal.
Regresión no lineal: Cuando f(X) no es una función lineal.
En tercer lugar, en función de la naturaleza de la relación que exista entre las dos variables:
La variable X puede ser la causa del valor de la variable Y.
Por ejemplo, en toxicología, si X = Dosis de la droga e Y = Mortalidad, la mortalidad se atribuye a la dosis administrada y no a otras causas.
Puede haber simplemente relación entre las dos variables.
Por ejemplo, en un estudio de medicina en que se estudian las variables X = Peso e Y = Altura de un grupo de individuos, puede haber relación entre las dos, aunque difícilmente una pueda considerarse causa de la otra.
Concepto de Regresión lineal y mínimos cuadrados
Los modelos lineales, relaccionan una salida con un conjunto de predictores de interás usando asunciones lineales. Son un subconjunto de los modelos de regresión, que son la herramienta más importante herramienta de análisis estadístico para un científico de datos. Un modelo de regresión lineal por tanto tendrá siempre la siguiente forma.
Y = A + BX + u
Y= variable dependiente o endógena
X= variable independiente o explicativa
A, B = parámetros fijos y desconocidos
Notación
- letras normales (i.e. \(X\), \(Y\)) = se usan para representar variables observadas
- La letras griegas (i.e. \(\mu\), \(\sigma\)) = se usan para representar variables desconocidas y que queremos estimar como la media o varianza de la población
- \(X_1, X_2, \ldots, X_n\) describe \(n\) puntos de datos
- \(\bar X\), \(\bar Y\) = medias observadas de las variables aleatorias \(X\) and \(Y\)
- \(\hat \beta_0\), \(\hat \beta_1\) = estimaciones de los valores reales \(\beta_0\) and \(\beta_1\)
Media de la muestra
- Media de la muestra se define como \[\bar X = \frac{1}{n}\sum_{i=1}^n X_i\]
- centrar la varible aleatoria es restar a todos sus valores la media de la muestra, por lo que tras el proceso la curva se centrará en el origen \[\tilde X_i = X_i - \bar X\]
- mean of \(\tilde X_i\) = 0
Desviación estándar y varianza
- varianza se define como \[ S^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar X)^2 = \frac{1}{n-1} \left( \sum_{i=1}^n X_i^2 - n \bar X ^ 2 \right) \Leftarrow \mbox{shortcut for calculation}\]
- desviación estándar se define como \(S = \sqrt{S^2}\)
- Al haber realizado la raiz cuadrada de la varianza, tenemos la desviación estándar en las mismas unidades que los datos
- escalar las variables aleatorias es dividir todos sus valores por la desviación estándar \(X_i / S\)
- standard deviation of \(X_i / S\) = 1
Normalizar
- Normalizar una variable aleatoria es centrarla y escalarla \[Z_i = \frac{X_i - \bar X}{s}\]
- empirical mean = 0, empirical standard deviation = 1
- Conseguimos que la distribución esté siempre centrada en 0 y que los datos muestren que están a unidades = # desviaciones estándar de la media original
- Ejemplo: \(Z_1 = 2\) significaría que punto Z está a 2 desviaciones estándar más allá de la media original
LA NORMALIZACIÓN PERMITE COMPARAR DATOS NO COMPARABLES INCLUSO QUE TENGAN DISTINTAS UNIDADES
Covarianza y Correlación
Tenemos \((X_i, Y_i)\) = como pareja de datos
- covarianza se define como \[
Cov(X, Y) = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar X) (Y_i - \bar Y) = \frac{1}{n-1}\left( \sum_{i=1}^n X_i Y_i - n \bar X \bar Y\right)
\]
- tiene como unidades \(X \veces\) unidades de \(Y\)
correlación se define como \[Cor(X, Y) = \frac{Cov(X, Y)}{S_x S_y}\] donde \(S_x\) y \(S_y\) son la estimación de las desviaciones estándar para las observaciones de \(X\) y de \(Y\) respectivamente.
El valor es efectivamente la covarianza estandarizada en una cantidad sin unidades
\(Cor(X, Y) = Cor(Y, X)\)
\(-1 \leq Cor(X, Y) \leq 1\)
\(Cor(X,Y) = 1\) and \(Cor(X, Y) = -1\) sólo cuando las observaciones \(X\) o \(Y\) caen perfectamente en una linea de pendiente positiva y negativa respectivamente
\(Cor(X, Y)\) mide la fortaleza de la relación lineal entre \(X\) e \(Y\), con una relacción mucho más fuerte cada vez que \(Cor(X,Y)\) se acerca a -1 o 1
\(Cor(X, Y) = 0\) implica que no hay relación lineal.
qownnotes-media-dXspNO
qownnotes-media-CwzJNE
qownnotes-media-OXitOI
Galton intentaba inferir la altura de los hijos a partir de la altura de los padres.
qownnotes-media-TlRIOL
{r galton,fig.height=3.5,fig.width=8}
library(UsingR); data(galton); library(reshape); long <- melt(galton)
g <- ggplot(long, aes(x = value, fill = variable))
g <- g + geom_histogram(colour = "black", binwidth=1)
g <- g + facet_grid(. ~ variable)
g
qownnotes-media-SqyrfP
Antes de ponernos a ver si podemos inferir la altura de los hijos por la altura de los padres a través del método de los mínimos cuadrados, explicaremos que el mínimo cuadrado es el valor mínimo de cada resta de un valor menos la media, todo ello al cuadrado.
qownnotes-media-ndPBIN
qownnotes-media-sjtBnE
# load necessary packages/install if needed
library(ggplot2); library(UsingR); data(galton)
# function to plot the histograms
myHist <- function(mu){
# calculate the mean squares
mse <- mean((galton$child - mu)^2)
# plot histogram
g <- ggplot(galton, aes(x = child)) + geom_histogram(fill = "salmon",
colour = "black", binwidth=1)
# add vertical line marking the center value mu
g <- g + geom_vline(xintercept = mu, size = 2)
g <- g + ggtitle(paste("mu = ", mu, ", MSE = ", round(mse, 2), sep = ""))
g
}
# manipulate allows the user to change the variable mu to see how the mean squares changes
# library(manipulate); manipulate(myHist(mu), mu = slider(62, 74, step = 0.5))]
# plot the correct graph
myHist(mean(galton$child))El resultado como podemos observar es que la media (lógicamente) es el valor que minimiza al máximo el valor del mínimo cuadrado, si elegimos otro valor como media el mse ya no será el mínimo.
Una vez sabemos el valor de mse, vamos a mostrar un gráfico de 928 puntos con las alturas de los padres y los hijos emparejados. Las alturas de las mujeres han sido ajustadas a 1.08 respecto a los hombre para poder hacer la media. En el plot usamos la función jitter para que no se superpongan mucho los valores coincidentes.
{r,dependson="galton",fig.height=4,fig.width=4, fig.align='center'}
ggplot(galton, aes(x = parent, y = jitter(child))) + geom_point(shape=1)
qownnotes-media-BUMWUG
Ahora lo que haremos será añadir una linea al gráfico de color azul. Esta línea representa la mínima variación de los datos alrededor de él. Si la pendiente de la recta es mayor que cero significa que la altura de los padres afecta a la altura de los hijos.
La pendiente de la recta es la estimación del coeficiente, o multiplicador de los padres, la variable independiente de nuestros datos. Con summary se puede ver la pendiente de la recta de regresión:
regrline <- lm(child ~ parent,galton)
summary(regrline)
Call: lm(formula = child ~ parent, data = galton)
Residuals: Min 1Q Median 3Q Max -7.8050 -1.3661 0.0487 1.6339 5.9264
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.94153 2.81088 8.517 <2e-16
parent 0.64629 0.04114 15.711 <2e-16
Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ‘’ 1
Residual standard error: 2.239 on 926 degrees of freedom Multiple R-squared: 0.2105, Adjusted R-squared: 0.2096 F-statistic: 246.8 on 1 and 926 DF, p-value: < 2.2e-16
Vamos a añadir ahora dos líneas azules para indicar las medias de las alturas de los padres y de los hijos . Cabe destacar que ambas líneas se cruzan en el mismo punto, justo donde está el punto medio de la recta.
La pendiente de la curva roja muestra cuánto cambia en la dirección vertical en un cambio de la dirección horizontal, vemos que padres que están 1 pulgada por encima de la media suelen tener hijos que sólo están a 0.65 pulgadas de la media, el triángulo verde indica este punto.
qownnotes-media-DBHmup
Concepto de Residuos
- Residuo, \(e_i\) = diferencia entre la salida observada y la que se ha estimado \[e_i = Y_i - \hat Y_i\]
- También, distancia vertical entre el dato observado y la linea de regresión
- Los mínimos cuadrados minimizan \(\sum_{i=1}^n e_i^2\)
- \(e_i\) puede ser interpretado como el estimador del error de regresión, \(\epsilon_i\)
- \(e_i\) también puede ser interpretado como la salida (\(Y\)) con asociación lineal del predictor (\(X\)) eliminada
- o, “Y ajustada a X”
- \(E[e_i] = 0\) \(\rightarrow\) esto es porque la media de los residuos se espera que sea 0 (se asume distribución gaussiana)
- Esta distribución de Gauss implica que el error NO está en correlación con ningún predictor
mean(fitModel$residuals)= devuelve la media de los residuos \(\rightarrow\) debe ser igual a 0cov(fit$residuals, predictors)= devuelve la covarianza (mide correlación) entre residuos y predictores \(\rightarrow\) debe ser igual a 0
- \(\sum_{i=1}^n e_i = 0\) (si el interceptor está incluido) and \(\sum_{i=1}^n e_i X_i = 0\) (si una variable regresora \(X_i\) está incluida)
- para modelos de regresión lineal estándar
- positive residuals = por encima de la linea
- negative residuals = por debajo
- Los residuos al ponerser en gráfica puden destacar una pobre convergencia del modelo
A continuación hay algunos ejemplos de cómo calcular los residuos, las distancias entre las alturas reales de los hijos y la estimada dada la línea de regresión. Como todas las líneas, está caracterizada por dos parámetros, una pendiente y un intercept, usaremos los mínimos cuadrados para proporcionar dos ecuaciones en dos datos desconocidos, de manera que podamos resolver estos dos parámetros
…
|====== | 6%
…
|========== | 9%
lm(child ~ parent,galton)
Call: lm(formula = child ~ parent, data = galton)
Coefficients: (Intercept) parent
23.9415 0.6463
fit <- lm(child ~ parent,galton)
|============= | 12%
summary(fit)
Call: lm(formula = child ~ parent, data = galton)
Residuals: Min 1Q Median 3Q Max -7.8050 -1.3661 0.0487 1.6339 5.9264
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.94153 2.81088 8.517 <2e-16 ***
parent 0.64629 0.04114 15.711 <2e-16
|========================== | 25%
…
|============================= | 28%
1: I haven’t the slightest idea. 2: mch = ic + slopemph 3: mph = ic + slopemch
Selection: 3
1: mph = ic + slopemch 2: I haven’t the slightest idea. 3: mch = ic + slopemph
Selection: 3
|================================ | 31%
…
|==================================== | 34%
…
|======================================= | 38%
…
|========================================== | 41%
…
|============================================== | 44%
…
|================================================= | 47%
ols.ic <- fit$coef[1]
|==================================================== | 50%
ols.slope <- fit$coef[2]
|======================================================= | 53%
…
|========================================================== | 56%
rhs -lhs [1] -1.264198e-09 -2.527486e-09 -3.801688e-09 1.261469e-09 2.522938e-09 3.767127e-09
lhs-rhs [1] 1.264198e-09 2.527486e-09 3.801688e-09 -1.261469e-09 -2.522938e-09 -3.767127e-09
|============================================================== | 59%
all.equal(lhs,rhs) [1] TRUE
|================================================================= | 62%
varchild <- var(galton$child)
varChild <- var(galton$child)
|==================================================================== | 66%
varRes <- var(fit$residuals)
|======================================================================== | 69%
varEst(est(slope = ols.slope,intercept = ols.ic)) Error: could not find function “varEst” varEST <- var(est(slope = ols.slope,intercept = ols.ic))
varEst <- var(est(ols.slope, ols.ic))
|=========================================================================== | 72%
all.equal(varChild,sum(varRes,varEst)) [1] TRUE
all.equal(varChild,varEst+varRes) [1] TRUE
|============================================================================== | 75%
…
|================================================================================= | 78%
1: is less than the variance of data 2: is unknown without actual data 3: is greater than the variance of data
Selection: 1
|==================================================================================== | 81%
…
|======================================================================================== | 84%
efit <- lm(accel ~ mag+dist, attenu)
|=========================================================================================== | 88%
mean(efit$residuals) [1] -1.785061e-18
|============================================================================================== | 91%
cov(efit\(residuals,efit\)mag) Error in cov(efit\(residuals, efit\)mag) : supply both ‘x’ and ‘y’ or a matrix-like ‘x’ cov(x = efit\(residuals,y = efit\)mag) Error in cov(x = efit\(residuals, y = efit\)mag) : supply both ‘x’ and ‘y’ or a matrix-like ‘x’ cov(x = efit\(residuals) Error in cov(x = efit\)residuals) : supply both ‘x’ and ‘y’ or a matrix-like ‘x’ cov() Error in is.data.frame(x) : argument “x” is missing, with no default efit$mag NULL
cov(efit\(residuals, attenu\)mag) [1] 5.338694e-17
|================================================================================================== | 94%
cov(efit\(residuals, attenu\)dist) [1] 5.253433e-16
Concepto de estimación por mínimos cuadrados
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:Hmisc':
##
## src, summarize
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Warning: Ignoring unknown aesthetics: show_guide
qownnotes-media-UYlNgi
qownnotes-media-gaybyM
qownnotes-media-BpebuK
qownnotes-media-kWwUpE
qownnotes-media-nBqHaI
qownnotes-media-CmQKYz
## Warning: Ignoring unknown aesthetics: show_guide
Como podemos ver, la linea de regresión resume la relación entre la altura de los padres (predictor) y la altura de los hijos (outcome)
Ejercicios
…
|================================= | 32%
qownnotes-media-yULohl
|================================================= | 47%
1: .64 2: .44 3: 5 4: .70
|======================================================= | 53%
1: .64 2: .66 3: 5.0 4: 44
Selection: 2
1: .64 2: 44 3: .66 4: 5.0
Selection: 4
|============================================================ | 58%
…
|================================================================== | 63%
cor(gch_nor,gpa_nor) [1] 0.4587624
|======================================================================= | 68%
1: It is the same. 2: It is bigger. 3: It is smaller.
Selection: 1
|============================================================================= | 74%
l_nor <- lm(gch_nor~gpa_nor)
|================================================================================== | 79%
1: I have no idea 2: 1. 3: The correlation of the 2 data sets
Selection: 3
|======================================================================================== | 84%
1: the same as the original 2: I have no idea 3: 1. 4: correlation(X,Y) * sd(X)/sd(Y)
Selection: 4
|============================================================================================= | 89%
qownnotes-media-LnkyEN
TEST 1
qownnotes-media-iTHEhs
qownnotes-media-bTsdUo
qownnotes-media-egzUYt
qownnotes-media-jSoMNj
qownnotes-media-rGCbis
qownnotes-media-CMiOzc
qownnotes-media-rPUzYm
qownnotes-media-UQjwkz
qownnotes-media-elvFIV
Aplicación de regresión lineal y mínimos cuadrados
Hasta ahora sólo hemos cosiderado hasta ahora la estimación, que es útil, pero también necesitamos saber cómo extender nuestras estimaciones a la población. Éste es el proceso de inferencia estadística. Nuestra aproximación a la inferencia estadística será a través de modelos estadísticos. Como mínimo, necesitamos algunas asunciones en la distribución de errores, nos focalizaremos en modelos sobre gausanidad.
Vamos a ver cómo hacer a través de mínimos cuadrados cálculos sobre al ecuación modificado el punto de corte con y (intercept) y la escala de la pendiente (por ejemplo porque tenga un valor muy alto).
Usamos el ejemplo del precio de los diamantes con respecto a sus quilates (carat)
Plot of the data
Intercept es el corte con X y el valor de carat es la pendiente de la recta, cada incremento de 1 kilate supone 3721$
fit <- lm(price ~ carat, data = diamond)
coef(fit)## (Intercept) carat
## -259.6259 3721.0249
qownnotes-media-jiuSWl
No queremos para nada saber el precio de 0 kilates (-259.6), por lo que lo que podemos calcular el precio de la media de kilates para el valor de intercept Para hacer un cálculo aritmétco dentro de una función hay que poner I( )
Getting a more interpretable intercept
fit2 <- lm(price ~ I(carat - mean(carat)), data = diamond)
coef(fit2)## (Intercept) I(carat - mean(carat))
## 500.0833 3721.0249
Ahora vemos que el cambio de 1 kilate es demasiado alto, por lo que nos puede interesar escalarlo para tener una mejor idea del valor
fit3 <- lm(price ~ I(carat * 10), data = diamond)
coef(fit3)## (Intercept) I(carat * 10)
## -259.6259 372.1025
qownnotes-media-FUyJBN
Nuestro fin último es poder estimar el valor de nuevos diamantes, en este caso queremos inferir el precio de 3 diamantes
Predicting the price of a diamond
newx <- c(0.16, 0.27, 0.34)
coef(fit)[1] + coef(fit)[2] * newx## [1] 335.7381 745.0508 1005.5225
predict(fit, newdata = data.frame(carat = newx))## 1 2 3
## 335.7381 745.0508 1005.5225
qownnotes-media-xCTQRV
Si hacemos predict(fit) nos devolverá únicamente los valores predecidos para las x del dataset
Podemos mostrar los valores que hemos predecido en la misma gráfica
Ejercicios
ones <- rep(1, nrow(galton))
lm(child ~ ones + parent -1, galton)
Call: lm(formula = child ~ ones + parent - 1, data = galton)
Coefficients: ones parent
23.9415 0.6463
lm(child ~ parent, galton)
Call: lm(formula = child ~ parent, data = galton)
Coefficients: (Intercept) parent
23.9415 0.6463
1: True 2: False
Selection: 1
qownnotes-media-bLHzWL
…
|======================================== | 35%
1: The variable itself 2: The outcome 3: The constant, 1
Selection: 3
|============================================= | 38%
lm(child ~ 1, galton)
Call: lm(formula = child ~ 1, data = galton)
Coefficients: (Intercept)
68.09
1: True 2: False
Selection: 1
head(trees) Constant Girth Height Volume 1 1 8.3 70 10.3 2 1 8.6 65 10.3 3 1 8.8 63 10.2 4 1 10.5 72 16.4 5 1 10.7 81 18.8 6 1 10.8 83 19.7
# Regress the given variable on the given predictor,
# suppressing the intercept, and return the residual.
regressOneOnOne <- function(predictor, other, dataframe){
# Point A. Create a formula such as Girth ~ Height -1
formula <- paste0(other, " ~ ", predictor, " - 1")
# Use the formula in a regression and return the residual.
resid(lm(formula, dataframe))
}
# Eliminate the specified predictor from the dataframe by
# regressing all other variables on that predictor
# and returning a data frame containing the residuals
# of those regressions.
eliminate <- function(predictor, dataframe){
# Find the names of all columns except the predictor.
others <- setdiff(names(dataframe), predictor)
# Calculate the residuals of each when regressed against the given predictor
temp <- sapply(others, function(other)regressOneOnOne(predictor, other, dataframe))
# sapply returns a matrix of residuals; convert to a data frame and return.
as.data.frame(temp)
}
1: Volume ~ Girth 2: Girth ~ Volume - 1 3: Volume ~ Girth - 1
Selection: 3
fit <- lm(Volume ~ Girth + Height + Constant -1, trees)
trees2 <- eliminate(“Girth”, trees)
head(trees2) Constant Height Volume 1 0.4057735 24.38809 -9.793826 2 0.3842954 17.73947 -10.520109 3 0.3699767 14.64038 -11.104298 4 0.2482677 14.29818 -9.019900 5 0.2339490 22.19910 -7.104089 6 0.2267896 23.64956 -6.446183
1: Computational precision was insufficient. 2: The constant, 1, has been replaced by its residual when regressed against Girth. 3: There must be some mistake
Selection: View(trees2) Enter an item from the menu, or 0 to exit Selection: 2
fit2 <- lm(Volume ~ Height + Constant-1, trees2)
lapply(list(fit, fit2), coef) [[1]] Girth Height Constant 4.7081605 0.3392512 -57.9876589
[[2]] Height Constant 0.3392512 -57.9876589
Call: lm(formula = Volume ~ Constant - 1, data = eliminate(“Height”, trees2))
Coefficients: Constant
-57.99
1: Subtract the mean from the outcome and each regressor. 2: Pick any regressor and replace the outcome and all other regressors by their residuals against the chosen one.
Selection: 2
Varianza de Residuos
Los residuos representan la variación que no ha podido ser explicada por nuestro modelo. Es importante notar la diferencia entre residuos y errores. Los errores son errores reales no observados desde los coeficientes conocidos, mientras que los residuos son errores observables desde los coeficientes observados. De alguna manera, los residuos son la estimación de los errores
Los residuos son la diferencia entre la estimación de las y y el valor real que en conjunto de datos. Tenemos que recordar que para que la regresión lineal sea aplicable, hay una serie de hipótesis como que la suma de los residuos es cero y que la varianza de los residuos es uniforme
Recordamos que si hacemos predict(fit) obtenemos los valores de la predicción de las x en lugar de su valor real, aquí verermos cómo la diferencia entre calcular los residuos con resid() y calculándolo a mano, la máxima diferencia entre ambos resultados es muy cercano a 0
qownnotes-media-pNWqxF
Además podemos comprobar una de nuestras hipótesis para poder aplcar la regresión lineal, y es que la suma de los reiduos es 0, así como si lo multiplicamos por x
qownnotes-media-kbAnjm
{r, echo = FALSE, fig.height=5, fig.width=5}
plot(diamond$carat, diamond$price,
xlab = "Mass (carats)",
ylab = "Price (SIN $)",
bg = "lightblue",
col = "black", cex = 2, pch = 21,frame = FALSE)
abline(fit, lwd = 2)
for (i in 1 : n)
lines(c(x[i], x[i]), c(y[i], yhat[i]), col = "red" , lwd = 2)
qownnotes-media-DuQGjF
Muchas veces los residuos nos permiten ver patrones en los datos, por lo que es muy útil ponerlo en el eje de las x
{r, echo = FALSE, fig.height=5, fig.width=5}
n =
plot(x, e,
xlab = "Mass (carats)",
ylab = "Residuals (SIN $)",
bg = "lightblue",
col = "black", cex = 2, pch = 21,frame = FALSE)
abline(h = 0, lwd = 2)
for (i in 1 : n)
lines(c(x[i], x[i]), c(e[i], 0), col = "red" , lwd = 2)
Este no es el caso ya que nuestro modelo se ajusta muy bien, pero por ejemplo nos ha servido para ver que hay muchos casos en los que para un mismo tamaño, no es siempre el mismo precio
qownnotes-media-QLKGdA
Hay veces, que aunque los datos no sean lineales, mostrar su recta de regresión es muy útil para detectar su patrón, aunque la recta o no nos sirva para hacer una predicción adecuada (the model doesnt fit)
No linear data (metemos sen() en la función y introduciendo también algo de ruido normalizado, en las x una distribución uniforme runif,)
{r, echo = FALSE, fig.height=5, fig.width=5}
x = runif(100, -3, 3); y = x + sin(x) + rnorm(100, sd = .2);
library(ggplot2)
g = ggplot(data.frame(x = x, y = y), aes(x = x, y = y))
g = g + geom_smooth(method = "lm", colour = "black")
g = g + geom_point(size = 7, colour = "black", alpha = 0.4)
g = g + geom_point(size = 5, colour = "red", alpha = 0.4)
g
qownnotes-media-NprLkt
Así a primera vista nos podría parecer un modelo acertado, pero si observamos de cerca sus residuos:
{r, echo = FALSE, fig.height=5, fig.width=5}
g = ggplot(data.frame(x = x, y = resid(lm(y ~ x))),
aes(x = x, y = y))
g = g + geom_hline(yintercept = 0, size = 2);
g = g + geom_point(size = 7, colour = "black", alpha = 0.4)
g = g + geom_point(size = 5, colour = "red", alpha = 0.4)
g = g + xlab("X") + ylab("Residual")
g
Vemos que realmente hay mucha variabilidad respecto a la recta de regresión
qownnotes-media-SuXmwZ
E incluso vemos el patrón mucho más claramente en los datos
qownnotes-media-KBfCUc
Otra de las cosas que podemos comprobar es si se cumple la segunda hipótesis, que es que la distribución de la varianza de los residuos sea uniforme. esto lo vamos a comprobar con la heterocedasticidad
qownnotes-media-uMHevD
Aquí aparentemente vemos que la recta es muy acertada
{r, echo = FALSE, fig.height=4.5, fig.width=4.5}
x <- runif(100, 0, 6); y <- x + rnorm(100, mean = 0, sd = .001 * x);
g = ggplot(data.frame(x = x, y = y), aes(x = x, y = y))
g = g + geom_smooth(method = "lm", colour = "black")
g = g + geom_point(size = 7, colour = "black", alpha = 0.4)
g = g + geom_point(size = 5, colour = "red", alpha = 0.4)
g
qownnotes-media-oYqzkP
Sin embargo si echamos un vistazo más profundo a los residuos, vemos que su distribución no es uniforme, si no que crece con las x, esto es lo que se llama heterocedasticidad, y lo conseguimos metiendo las x dentro de la función rnorm de las y
y <- x + rnorm(100, mean = 0, sd = .001 * x)
{r, echo = FALSE, fig.height=4.5, fig.width=4.5}
g = ggplot(data.frame(x = x, y = resid(lm(y ~ x))),
aes(x = x, y = y))
g = g + geom_hline(yintercept = 0, size = 2);
g = g + geom_point(size = 7, colour = "black", alpha = 0.4)
g = g + geom_point(size = 5, colour = "red", alpha = 0.4)
g = g + xlab("X") + ylab("Residual")
g
qownnotes-media-NyOltz
Volviendo al ejempo de los diamantes, vamos a añadir al dataset los residuos de cada par de valores, y vamos a dibujarlos frente a los kilates, los residuos tendrán siempre las unidades de las y
{r, echo = FALSE, fig.height=4.5, fig.width=4.5}
diamond$e <- resid(lm(price ~ carat, data = diamond))
g = ggplot(diamond, aes(x = carat, y = e))
g = g + xlab("Mass (carats)")
g = g + ylab("Residual price (SIN $)")
g = g + geom_hline(yintercept = 0, size = 2)
g = g + geom_point(size = 7, colour = "black", alpha=0.5)
g = g + geom_point(size = 5, colour = "blue", alpha=0.2)
g
qownnotes-media-ByaBpa
Lo primero que hacemos es calcular los residuos frente a la media del precio (esto es frente al interception), el otro vector es ls residuos frente a la recta de regresión
Luego creo dos facores, uno de interception y otro de inter y slope
{r, echo = FALSE, fig.height=4.5, fig.width=4.5}
e = c(resid(lm(price ~ 1, data = diamond)),
resid(lm(price ~ carat, data = diamond)))
fit = factor(c(rep("Itc", nrow(diamond)),
rep("Itc, slope", nrow(diamond))))
g = ggplot(data.frame(e = e, fit = fit), aes(y = e, x = fit, fill = fit))
g = g + geom_dotplot(binaxis = "y", size = 2, stackdir = "center", binwidth = 20)
g = g + xlab("Fitting approach")
g = g + ylab("Residual price")
g
La varianza de los residuos por tanto será la suma al cuadrado de los residuos dividido n
qownnotes-media-VHpgFq
La razón por la que se resta 2 es por quitarle 2 grados de libertad, uno por conocer el final (intercerp) como con la varianza y otro grado de libertad por la covarianza entre x e y
qownnotes-media-QdNeXf
En el primer gráfico vemos la variación en el precio de los diamantes alrededor de la media del precio.
En el gráfico azul vemos la variación frente a la linea de regresión
qownnotes-media-xBrfKR
qownnotes-media-lbWthX
qownnotes-media-fWNVhO
R^2
Es el porcentaje total de la variabilidad explicada por la relacción lineal con el predictor
qownnotes-media-ALluas
En los siguientes ejemplos vamos a ver cómo a pesar de que todo indica una buena realacción lineal, en realidad no la hay
qownnotes-media-ISIpjA
En el primer gráfico se ve que hay mucho ruido, en el segundo que es una curva que no se ajusta, en el tercero hay un outlier que interfiere mucho y en el cuarto hay un punto que hace que se degenere.
qownnotes-media-hvZcHw
Ejercicios
lm(child ~ parent,galton)
Call: lm(formula = child ~ parent, data = galton)
Coefficients: (Intercept) parent
23.9415 0.6463
fit <- lm(child ~ parent, galton)
sqrt(sum(fit$residuals^2) / (n - 2)) [1] 2.238547
summary(fit)$sigma [1] 2.238547
sqrt(deviance(fit)/(n-2)) [1] 2.238547
1: Yi_hat-mean(Yi) 2: Yi-mean(Yi) 3: Yi-Yi_hat
Selection: 2
1: Yi_hat-mean(Yi) 2: Yi-Yi_hat 3: Yi-mean(Yi)
Selection: 2
mu <- mean(galton$child)
sTot <- sum((galton$child-mu)^2)
sRes <- deviance(fit)
1 - sRes/sTot [1] 0.2104629
summary(fit)$r.squared [1] 0.2104629
cor(galton\(child,galton\)parent)^2 [1] 0.2104629
Inferencia
Es el proceso de sacar conclusiones acerca de la población desde una muestra. En inferencia estadística, debemos ser siempre conscientes de la incertidumbre en nuestras estimaciones concienzudamente. Los test de hipótesis e intervalos de confianza son las formas más comunes de inferencia estadística.
qownnotes-media-SZlQQz
Podemos calcular todo a mano para ver que todo cuadra
{r}
library(UsingR); data(diamond)
y <- diamond$price; x <- diamond$carat; n <- length(y)
beta1 <- cor(y, x) * sd(y) / sd(x)
beta0 <- mean(y) - beta1 * mean(x)
e <- y - beta0 - beta1 * x
sigma <- sqrt(sum(e^2) / (n-2))
ssx <- sum((x - mean(x))^2)
seBeta0 <- (1 / n + mean(x) ^ 2 / ssx) ^ .5 * sigma
seBeta1 <- sigma / sqrt(ssx)
tBeta0 <- beta0 / seBeta0; tBeta1 <- beta1 / seBeta1
pBeta0 <- 2 * pt(abs(tBeta0), df = n - 2, lower.tail = FALSE)
pBeta1 <- 2 * pt(abs(tBeta1), df = n - 2, lower.tail = FALSE)
coefTable <- rbind(c(beta0, seBeta0, tBeta0, pBeta0), c(beta1, seBeta1, tBeta1, pBeta1))
colnames(coefTable) <- c("Estimate", "Std. Error", "t value", "P(>|t|)")
rownames(coefTable) <- c("(Intercept)", "x")
Comparamos los resultados del calculo a mano y el calculo con lm
{r}
coefTable
fit <- lm(y ~ x);
summary(fit)$coefficients
qownnotes-media-PhOFwb
Y finalmente hacemos el cálculo del intervalo de confianza de nuestra estimación
qownnotes-media-bsQtRI
Podemos decir por tanto que con un 95% de probabilidad, el incremento de precio sobre 0.1 kilates en diamantes estará entre 355.63 y 388.56
Vamos a generar predicciones para nuevos valores de X al azar
{r, fig.height=5, fig.width==5, echo = FALSE, results='hide'}
library(ggplot2)
newx = data.frame(x = seq(min(x), max(x), length = 100))
p1 = data.frame(predict(fit, newdata= newx,interval = ("confidence")))
p2 = data.frame(predict(fit, newdata = newx,interval = ("prediction")))
p1$interval = "confidence"
p2$interval = "prediction"
p1$x = newx$x
p2$x = newx$x
dat = rbind(p1, p2)
names(dat)[1] = "y"
g = ggplot(dat, aes(x = x, y = y))
g = g + geom_ribbon(aes(ymin = lwr, ymax = upr, fill = interval), alpha = 0.2)
g = g + geom_line()
g = g + geom_point(data = data.frame(x = x, y=y), aes(x = x, y = y), size = 4)
g
qownnotes-media-HxdQrL
qownnotes-media-zSGDCG
all <- lm(Fertility ~ ., swiss)
summary(all)
Call: lm(formula = Fertility ~ ., data = swiss)
Residuals: Min 1Q Median 3Q Max -15.2743 -5.2617 0.5032 4.1198 15.3213
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 66.91518 10.70604 6.250 1.91e-07 Agriculture -0.17211 0.07030 -2.448 0.01873
Examination -0.25801 0.25388 -1.016 0.31546
Education -0.87094 0.18303 -4.758 2.43e-05 * Catholic 0.10412 0.03526 2.953 0.00519 Infant.Mortality 1.07705 0.38172 2.822 0.00734
Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ‘’ 1
Residual standard error: 7.165 on 41 degrees of freedom Multiple R-squared: 0.7067, Adjusted R-squared: 0.671 F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10
1: 0.1 2: 0.05 3: 0.01 4: R doesn’t say
Selection: Enter an item from the menu, or 0 to exit Selection: 2
all <- lm(Fertility ~ Agriculture, swiss)
summary(lm(Fertility ~ Agriculture, swiss))
Call: lm(formula = Fertility ~ Agriculture, data = swiss)
Residuals: Min 1Q Median 3Q Max -25.5374 -7.8685 -0.6362 9.0464 24.4858
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 60.30438 4.25126 14.185 <2e-16 ** Agriculture 0.19420 0.07671 2.532 0.0149
Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ‘’ 1
Residual standard error: 11.82 on 45 degrees of freedom Multiple R-squared: 0.1247, Adjusted R-squared: 0.1052 F-statistic: 6.409 on 1 and 45 DF, p-value: 0.01492
1: 0.19420 2: 60.30438 3: 0.07671 4: *
Selection: 1
1: They would be uncorrelated 2: They would be correlated 3: I would not be able to guess without more information
Selection: 2
cor(swiss\(Examination,swiss\)Education) [1] 0.6984153
cor(swiss\(Agriculture,swiss\)Education) [1] -0.6395225
makelms <- function(){
# Store the coefficient of linear models with different independent variables
cf <- c(coef(lm(Fertility ~ Agriculture, swiss))[2],
coef(lm(Fertility ~ Agriculture + Catholic,swiss))[2],
coef(lm(Fertility ~ Agriculture + Catholic + Education,swiss))[2],
coef(lm(Fertility ~ Agriculture + Catholic + Education + Examination,swiss))[2],
coef(lm(Fertility ~ Agriculture + Catholic + Education + Examination +Infant.Mortality, swiss))[2])
print(cf)
}
# Regressor generation process 1.
rgp1 <- function(){
print("Processing. Please wait.")
# number of samples per simulation
n <- 100
# number of simulations
nosim <- 1000
# set seed for reproducability
set.seed(4321)
# Point A:
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
# Point B:
betas <- sapply(1 : nosim, function(i)makelms(x1, x2, x3))
round(apply(betas, 1, var), 5)
}
# Regressor generation process 2.
rgp2 <- function(){
print("Processing. Please wait.")
# number of samples per simulation
n <- 100
# number of simulations
nosim <- 1000
# set seed for reproducability
set.seed(4321)
# Point C:
x1 <- rnorm(n)
x2 <- x1/sqrt(2) + rnorm(n) /sqrt(2)
x3 <- x1 * 0.95 + rnorm(n) * sqrt(1 - 0.95^2)
# Point D:
betas <- sapply(1 : nosim, function(i)makelms(x1, x2, x3))
round(apply(betas, 1, var), 5)
}
TEST2
qownnotes-media-FATtcx
qownnotes-media-FRaEUb
qownnotes-media-vvKXCO
qownnotes-media-vvjBvz
qownnotes-media-ytKPiJ
qownnotes-media-GXWwfn
qownnotes-media-BwKEvy
qownnotes-media-NZjjlP
Regresión multivariable
Cuando se hace un análisis multivariables, es muy importante escoger qué variables van a intervenir en el modelo, y de esa manera eliminar variables de confusión que pueden hacer que la regresón que estamos llevando a cabo sea errónea.
Vamos a comparar una regresión entre la fertilidad comparada con el resto de variables, nos vamos a fijar en su relacción con la agricultura, y luego vamos a ver su relacción excluyendo el resto de variables para ver cómo difiere. Esta todo en porcentajes.
qownnotes-media-AOeeST
qownnotes-media-uZewce
Se puede hacer gráficamente de dos maneras
qownnotes-media-KWvukl
qownnotes-media-NjyqIW
qownnotes-media-dIhRXq
Vamos a hacerlo con lm
qownnotes-media-AjSkFc
qownnotes-media-FXaRjc
El t value no es más que dividir el estimate/std Errror, y la última columna nos da un p value two sided que nos indica si es estadísticamente significativo
Ahora vamos a analizar la fertilidad y agricultura exclusivamente
qownnotes-media-rnBQSz
Los resultados son muy distintos, pero incluso me sigue diciendo que es estadisticamente significativo.
Este estudio es algo dinámico, tendremos que ir quitando y poniendo varables para probar su efecto en el modelo.
qownnotes-media-nplXCb
lm normalmente nos devuelva NA es porque esa variable es una combinación lineal de otras y no aporta ningún valor al modelo.
COMPARAR AGRUPACIONES EN VARIABLE CATEGORICA
Ahora vamos a ver cómo se puede usar lm para comparar resultados frente a una variable factor
En este caso vamos a usar el dataset de insectos, con el spray correspondiente y el número de insectors que mata
qownnotes-media-OvnTJV
qownnotes-media-UBmJtf
Si no indicamos nada, se toma como referencia del factor el primer grupo, en este caso A, que actuará como Intercep
Linear model fit, group A is the reference
{r, echo= TRUE}
summary(lm(count ~ spray, data = InsectSprays))$coef
Aquí lo que estamos obteniendo son el el coficiente la diferencia con la media referencia (intercept) y el Pr(t) nos dirá si la diferencia entre los grupos es estadísticamente significativa (vemos en el caso de A y B que no, por ejemplo, como se puede ver en el boxplot).
Aquí estamos asumiendo que estamos tratando con una distribución normal y que la varianza entre los grupos es homogénea, cosa que no es cierta (los valores son siempre mayores que 0) por lo que aquí el pvalue de ttest no será válido, pero nos sirve para ver el ejemplo. Para contar datos será más apropiada Poisson GLM que veremos más adelante.
qownnotes-media-ffKRON
Siempre tenemos que comparar con algo, si queremos que sea otra la referencia tenemos que hacer relevel
Reordering the levels
{r}
spray2 <- relevel(InsectSprays$spray, "C")
summary(lm(count ~ spray2, data = InsectSprays))$coef
qownnotes-media-DpIsru
Y qué pasaría si omitimos el intercept? En este caso no compararíamos los grupos entre sí y lo estaríamos comparando con 0 como referencia, por tanto los pvalues sería si la diferencia con 0 es estadísticamente significativa, y en los coeficientes únicamente las medias de cada grupo.
What if we omit the intercept?
{r, echo= TRUE}
summary(lm(count ~ spray - 1, data = InsectSprays))$coef
unique(ave(InsectSprays$count, InsectSprays$spray))
qownnotes-media-xINqMe
Líneas de regresión por grupos
Vamos a investigar qué pasa si dos variables que estamos usando interfieren una con la otra. Para simplificarlo usaremos una variable Bimodal en base a la religion, 1 catolico 0 protestante, y vamos a sacar una línea de regresión por cada tipo.
qownnotes-media-KcMfBL
qownnotes-media-kSzOJt
Vamos a sacar una linea de regresión para las provincias con mayoría católica y otra línea para las provincias con mayoría protestante.
Si las variables no interfieren entre ellas tendremos 2 líneas con distinto intercept y misma pendiente:
Modelo 1 .- no se tiene en cuenta religion
qownnotes-media-UuDtzx
qownnotes-media-yxLdSz
qownnotes-media-jSVfeC
Modelo 2 .- Los grupos de religión no interfieren entre sí
qownnotes-media-LojhDs
qownnotes-media-RsIVbn
Especificamos que sea factor, porque hay veces que los 0,1,2 los identifica como coeficientes crecientes, que cada valor corresponde al doble del inmediatamente anterior
qownnotes-media-sHVqoU
qownnotes-media-CYrSAu
Modelo 3.- Los grupos de religión sí interfieren
qownnotes-media-lnCMeq
qownnotes-media-VKTKln
Esto lo hacemos con el * , de modo que el primer intercept será de los no catolicos, y el tercer estimador será el intercept de los catolicos, mientras que el segundo y el cuarto serán las pendientes de los no catolicos y catolicos
qownnotes-media-XFrChM
qownnotes-media-nokyRe
…
|======= | 5%
dim(hunger) [1] 948 13
|=========== | 8%
948 [1] 948
|============== | 11%
names(hunger) [1] “X” “Indicator” “Data.Source” “PUBLISH.STATES” “Year” “WHO.region” “Country” “Sex”
[9] “Display.Value” “Numeric” “Low” “High” “Comments”
|================== | 14%
…
|===================== | 16%
fit <- lm(Numeric ~ Year,hunger)
|========================= | 19%
summary(fit)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 634.479660 121.1445995 5.237375 2.007699e-07 Year -0.308397 0.0605292 -5.095012 4.209412e-07
|============================ | 22%
1: -0.30840 2: 634.47966 3: 121.14460 4: 0.06053
Selection: 1
|================================ | 24%
1: As time goes on, the rate of hunger increases 2: I haven’t a clue 3: As time goes on, the rate of hunger decreases
Selection: 3
|=================================== | 27%
1: the percentage of hungry children at year 0 2: the number of hungry children at year 0 3: the number of children questioned in the survey
Selection: 1
|======================================= | 30%
fit <- lm(Numeric ~ Year,hunger[hunger$Sex==“Female”]) Error in
[.data.frame(hunger, hunger\(Sex == "Female") : undefined columns selected fit <- lm(Numeric ~ Year,hunger[,hunger\)Sex==“Female”]) Error in[.data.frame(hunger, , hunger$Sex == “Female”) : undefined columns selected fit <- lm(Numeric ~ Year,hunger)
lmF <- lm(hunger\(Numeric[hunger\)Sex==“Female”] ~ hunger\(Year[hunger\)Sex==“Female”])
|========================================== | 32%
lmM <- lm(hunger\(Numeric[hunger\)Sex==“Male”] ~ hunger\(Year[hunger\)Sex==“Male”])
|============================================== | 35%
qownnotes-media-tBfMNp
1: lmM\(coef[1] 2: lmF\)coef[2] 3: lmM$coef[2]
Selection: 3
|===================================================== | 41%
…
|======================================================== | 43%
lmBoth <- lm(Numeric ~ Year+Sex,hunger)
|============================================================ | 46%
summary(lmBoth)
Call:
lm(formula = Numeric ~ Year + Sex, data = hunger)
Residuals:
Min 1Q Median 3Q Max
-25.472 -11.297 -1.848 7.058 45.990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 633.5283 120.8950 5.240 1.98e-07 ***
Year -0.3084 0.0604 -5.106 3.99e-07 ***
SexMale 1.9027 0.8576 2.219 0.0267 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.2 on 945 degrees of freedom
Multiple R-squared: 0.03175, Adjusted R-squared: 0.0297
F-statistic: 15.49 on 2 and 945 DF, p-value: 2.392e-07
|=============================================================== | 49%
…
|=================================================================== | 51%
1: 635.431 2: 633.2199 3: 1.9027 4: I can’t tell since the data starts at 1970.
Selection: 3
1: I can’t tell since the data starts at 1970. 2: 635.431 3: 633.2199 4: 1.9027
Selection: 2
|====================================================================== | 54%
1: the annual decrease in percentage of hungry males 2: the annual decrease in percentage of hungry children of both genders 3: the annual decrease in percentage of hungry females
Selection: 2
|========================================================================== | 57%
qownnotes-media-BCnmqs
|============================================================================= | 59%
1: they have the same slope 2: I have no idea 3: they have slopes that are very close
Selection: 1
|================================================================================= | 62%
…
|==================================================================================== | 65%
lmInter <- lm(Numeric ~ Year,Sex,Sex*Year,hunger) Error in is.data.frame(data) : object ‘Sex’ not found lmBoth <- lm(Numeric ~ Year+Sex,hunger)
lmInter <- lm(Numeric ~ Year + Sex +Year*Sex, hunger)
|======================================================================================== | 68%
summary(lmInter)
Call:
lm(formula = Numeric ~ Year + Sex + Year * Sex, data = hunger)
Residuals:
Min 1Q Median 3Q Max
-25.913 -11.248 -1.853 7.087 46.146
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 603.50580 171.05519 3.528 0.000439 ***
Year -0.29340 0.08547 -3.433 0.000623 ***
SexMale 61.94772 241.90858 0.256 0.797946
Year:SexMale -0.03000 0.12087 -0.248 0.804022
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.21 on 944 degrees of freedom
Multiple R-squared: 0.03181, Adjusted R-squared: 0.02874
F-statistic: 10.34 on 3 and 944 DF, p-value: 1.064e-06
|=========================================================================================== | 70%
1: The model doesn’t say. 2: 61.94772 3: 603.5058
Selection: 3
|=============================================================================================== | 73%
1: 61.94772 2: 665.4535 3: 603.5058 4: The model doesn’t say.
Selection: 2
|================================================================================================== | 76%
1: -0.03000 2: 0.08547 3: The model doesn’t say. 4: -0.29340
Selection: 4
|====================================================================================================== | 78%
1: 0.12087 2: -0.03000 3: -0.32340 4: The model doesn’t say.
Selection: 2
1: The model doesn’t say. 2: 0.12087 3: -0.03000 4: -0.32340
Selection: 4
|========================================================================================================= | 81%
qownnotes-media-yXmRTO
|============================================================================================================= | 84%
1: Female 2: Male 3: They look about the same
Selection: 1
1: They look about the same 2: Male 3: Female
Selection: 2
|================================================================================================================ | 86%
…
|==================================================================================================================== | 89%
…
|======================================================================================================================= | 92%
…
|=========================================================================================================================== | 95%
1: b0+b2 2: b2+b3Y 3: b1+5b3 4: b2+b3*5
Selection: 4
|============================================================================================================================== | 97%
Ajustes y relacción entre variables
El Ajuste es la idea de poner regresores dentro del modelo lineal e investigar el role de una tercera variable en la relación entre los otros dos, ya que es habitual el caso en el que esta tercera variable pueda causar distorsión o confusión en las otras
Como ejemplo, podemos considerar los porcentajes de cancer de pulmón y los usos de pastillas de menta. Estamos mirando el volumen expiratorio forzado y el consumo de pastillas de menta. Si encontráramos significancia estadística en la relacción de la regresión, no sería inteligente ir rápidamente a los periódicos y decir que el consumo de pastillas de menta baja la capacidad pulmonar. En primer lugar, auqnue la relacción es sólida, no se sabe si es causal. Pero más importante en este caso, la verdadera razón son los hábitos tabáquicos. Fumar está muy relaccionado con ambas cosas, capacidad pulmonar y consumo de pastillas de menta. Cómo te defiendes de la acusación de que es variabilidad en los hábitos como fumador?
Si el hallazgo apareciera entre fumadores y no fumadores de manera separada, entonces se podría tener algo. La gente no lo creería a menos que eliminásemos la constante fumador. Esta es la idea de añadir una variable de regresión como ajuste. El coeficiente de interés es interpretado como el efecto del predictor en la respuesta, manteniendo ajuste de la variable constante.
Buscar los ajustes es una tarea complicada, ya que se trata de ir haciendo pruebas y comprobando resultados, muchas veces teniendo que recurrir a un experto en el sector para que nos aconseje.
A continuación camos a ver mediante ejemplos, cómo puede influir o no un ajuste mediante la comparación de dos grupos, digamos que unos reciben tratamiento y otros no lo reciben.
EJEMPLO 1
Las líneas horizontales son las medias de los grupos, mientras que el resto son las regresiones lineales teniendo en cuenta los grupos y sin tenerlos en cuenta.
qownnotes-media-DJgAwG
EJEMPLO2
qownnotes-media-CPMPro
EJEMPLO 3
La asociación marginal es la diferencia entre las medias
qownnotes-media-mkldqL
EJEMPLO 4
qownnotes-media-Tjioqp
EJEMPLO 5
qownnotes-media-uRllAi
EJEMPLO 6
Metemos una tercera variable representada por el color
qownnotes-media-GwKpZj
Si representáramos estos puntos en 3 dimensiones, nos daríamos cuenta que la Z forma un plano perfectamente alineado
Residuos y Diagnósticos de outliers
Como hemos visto con anterioridad, los residuos son las diferencias entre los valores del dataset y los calculados por nuestro modelo de regresión. Estudiar estos residuos nos va a ayudar a ver si nuestro modelo es bueno o no
Para explicar la influencia de determinados puntos en el modelo vamos a exponer distintas situaciones en las que tendríamos que decidir si un punto es un outlier y sacarlo del modelo o no. Partimos de una serie de puntos con una clara relacción lineal y una serie de puntos rojos sobre los que vamos a ver cómo impactaría su inclusión en el dataset. Vamos a analizar el peso y la capacidad que tiene para cambiar el modelo con ese peso (influencia). Se puede pensar en esto como si fuera un balancín, el mismo peso tendrá un efecto distinto si se posa cerca del eje que si lo ponemos en el extremo del balancín.
qownnotes-media-jCNAIa
Abajo izquierda: Está justo en medio de todos los puntos por lo que tiene nulo peso al estar justo en el centro y va a tener nula influencia al estar perfectamente alineado
Arriba a la derecha: Tiene mucho peso por al estar tan alejado del centro, pero al estar perfectamente aliineado no tendá apenas influencia.
Abajo derecha: Tendría mucha influencia y mucho peso , ya que está muy alejado del centro y no está alineado en absoluto.
Arriba izquierda: En este caso está totalmente alineado con el centro pero muy desalineado con la recta de regresión, con lo que no tendrá tanta influencia.
qownnotes-media-ZaMwrE
qownnotes-media-CSVMBA
R propociona diversas herramientas para hacer este análisis:
Lo malo de los residuos es que tienen unidad de medida, la misma que las Y, así que las dos primeras herramientas lo que hacer es estandarizarlas.
qownnotes-media-PuPqrM
qownnotes-media-skzXBX
Vamos a estudiar casos:
Caso 1
qownnotes-media-TyISYc
qownnotes-media-anTkny
qownnotes-media-VaVlKp
Vemos que en este caso df betas del punto 10,10 es mucho más largo que el resto al estar muy poco alineado con el resto, así como el hatvalue que nos indica que esta muy alejado del centro, con lo que aquí decidiríamos que es un outlier. Tendría mucha influencia
Caso 2
qownnotes-media-PYgHwp
qownnotes-media-ypsYyK
Aquí vemos que las dfbetas no es muy largo con respecto al resto al estar muy bien alineado, sin embargo el hatvalue sí que es muy grande, al estar muy lejos del centro. No tendá demasiada influencia
Existe una funcionalidad en R que nos ayuda a visualizar varias cosas acerca de los residuos y que nos ayudan a decidir si el modelo encaja o no es acertado:
qownnotes-media-oRbgEj
El primer cuadro: Aquí vemos si hay patrones en los residuos respecto los valores estimados, podríamos ver por ejemplo heteroscaticidad
The simplest diagnostic plot displays residuals versus fitted values. Residuals should be uncorrelated with the | fit, independent and (almost) identically distributed with mean zero. In this case there is a linear pattern involving all but one residual and the fit.
plot(lm(y ~ x, out2), which=1)
qownnotes-media-UjuEUC
Our influential outlier is in row 1 of the data. To exclude it is just a matter using out2[-1, ] rather than out2 as data
plot(lm(y ~ x, out2[-1,]), which=1)
qownnotes-media-QuRHSf
This plot has none of the patterned appearance of the first. It looks as we would expect if residuals were independently and (almost) identically distributed with zero mean, and were uncorrelated with the fit.
The change which inclusion or exclusion of a sample induces in coefficents is a simple measure of its influence. Subtract coef(fitno) from coef(fit) to see the change induced by including the influential first sample.
coef((lm(y ~ x, out2)) - coef((lm(y ~ x, out2[-1,])) (Intercept) x -0.01167866 -0.53363019
qownnotes-media-sSRaTe
qownnotes-media-PJEUic
qownnotes-media-nhIAkH
resno <- out2[1, “y”] - predict(fitno, out2[1,])
1-resid(fit)[1]/resno 1 0.6311547
qownnotes-media-QbzaEE
qownnotes-media-lkikYx
sigma <- sqrt(deviance(fit)/df.residual(fit))
qownnotes-media-wNuUSr
rstd <- resid(fit)/(sigma * sqrt(1-hatvalues(fit)))
qownnotes-media-UrFVcY
El segundo cuadro: Se usa para evaluar si los errorres siguen una distribución normal
qownnotes-media-bpEAwS
qownnotes-media-DIQlWY
Look at the outlier’s standardized residual, labeled 1 on the Normal QQ plot, it is almost 5 standar deviation from the mean
El tercer cuadro: Es igual que el primero pero con los datos escalados, lo que nos permitiría compararlo con otros directamente in preocuparnos de las undades
qownnotes-media-BVZMqo
qownnotes-media-bxkQXo
qownnotes-media-uAIqZk
El cuarto cuadro: Este nos muestra el peso frente a la influencia de cada punto
sigma1 <- sqrt(deviance(fitno)/df.residual(fitno))
Calculate the Studentized residual for the outlier sample by dividing resid(fit)[1] by the product of sigma1 and | sqrt(1-hatvalues(fit)[1]). There is no need to store this in a variable.
qownnotes-media-OjOrwb
qownnotes-media-cDYWdG
qownnotes-media-oQgkYI
qownnotes-media-sNOUMr
Selección del Modelo (VIF y ANOVA)
Para elegir el modelo correcto de predicción y elegir las variables adecuadas no existe una regla, habrá que ir mirado cada caso e ir tomando decisiones.
Quitar varables importantes, nos va dar resultados sesgados, mientras que meter variables de más no va a inferir en los resultados de la predicción pero lo que produce es incrementar el error estándar.
qownnotes-media-NkyPQS
En este primer caso vemos que añadir variables al modelo, no impacta demasido en los resultados, nos vamos a basar en el primer método que es VIF
VIF
VARIANCE INFLATION FACTORS
Es el incremento de la varianza por el regresor i comparado con la situación ideal donde es ortogonal con el resto de regresores. La raíz cuadrada del VIF es el incremento en el sd. Hay que recordar la inflación de la varianza es sólo una parte del cuadro, queremos incluir ciertas variables incluso si aumentan dramáticamente nuestra varianza.
qownnotes-media-iQcjNt
En este segundo ejempo en el que ponemos una gran dependencia entre la variables se ve claramente cómo varia:
qownnotes-media-KsHtiM
qownnotes-media-kkKeGA
qownnotes-media-mmNOIM
qownnotes-media-ZOcqNv
En cuanto al aumento del error estándar:
qownnotes-media-fXblbs
Ejercicios
In modeling, our interest lies in parsimonious, interpretable representations of the | data that enhance our understanding of the phenomena under study. Omitting variables | results in bias in the coefficients of interest - unless their regressors are | uncorrelated with the omitted ones. On the other hand, including any new variables | increases (actual, not estimated) standard errors of other regressors. So we don’t | want to idly throw variables into the model. This lesson is about the second of | these two issues, which is known as variance inflation.
makelms <- function(x1, x2, x3){
# Simulate a dependent variable, y, as x1
# plus a normally distributed error of mean 0 and
# standard deviation .3.
y <- x1 + rnorm(length(x1), sd = .3)
# Find the coefficient of x1 in 3 nested linear
# models, the first including only the predictor x1,
# the second x1 and x2, the third x1, x2, and x3.
c(coef(lm(y ~ x1))[2],
coef(lm(y ~ x1 + x2))[2],
coef(lm(y ~ x1 + x2 + x3))[2])
}
1: The coefficient of x2. 2: The coefficient of the intercept. 3: The coefficient of x1.
Selection: 3
1: x1 and x2 2: x1 3: x1, x2, and x3
Selection: 2
|============ | 21%
# Regressor generation process 1.
rgp1 <- function(){
print("Processing. Please wait.")
# number of samples per simulation
n <- 100
# number of simulations
nosim <- 1000
# set seed for reproducibility
set.seed(4321)
# Point A
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
# Point B
betas <- sapply(1 : nosim, function(i)makelms(x1, x2, x3))
round(apply(betas, 1, var), 5)
}
# Regressor generation process 2.
rgp2 <- function(){
print("Processing. Please wait.")
# number of samples per simulation
n <- 100
# number of simulations
nosim <- 1000
# set seed for reproducibility
set.seed(4321)
# Point C
x1 <- rnorm(n)
x2 <- x1/sqrt(2) + rnorm(n) /sqrt(2)
x3 <- x1 * 0.95 + rnorm(n) * sqrt(1 - 0.95^2)
# Point D
betas <- sapply(1 : nosim, function(i)makelms(x1, x2, x3))
round(apply(betas, 1, var), 5)
}
1: x1, x2, and x3 are correlated in both rgp1() and rgp2(). 2: x1, x2, and x3 are uncorrelated in both rgp1() and rgp2(). 3: x1, x2, and x3 are uncorrelated in rgp1(), but not in rgp2(). 4: x1, x2, and x3 are correlated in rgp1(), but not in rgp2().
Selection: 2
|============== | 25%
1: Computes the variance of each row. 2: Computes the variance of each column.
Selection: 1
rgp1() [1] “Processing. Please wait.” x1 x1 x1 0.00110 0.00111 0.00112
1: x2 2: x3
Selection: 2
|===================== | 38%
rgp2() [1] “Processing. Please wait.” x1 x1 x1 0.00110 0.00240 0.00981
In this case, variance inflation due to correlated regressors is clear, and is most | pronounced in the third model, y ~ x1 + x2 + x3, since x3 is the regressor most | strongly correlated with x1.
…
|========================== | 46%
…
|=============================== | 54%
qownnotes-media-IhrIMw
|================================= | 58%
mdl <- lm(fertility ~ Agriculture+Examination+Education+Catholic+Infant.Mortality,swiss) Error in eval(expr, envir, enclos) : object ‘fertility’ not found mdl <- lm(Fertility ~ Agriculture+Examination+Education+Catholic+Infant.Mortality,swiss)
|==================================== | 62%
vif(mdl) Agriculture Examination Education Catholic Infant.Mortality 2.284129 3.675420 2.774943 1.937160 1.107542
|====================================== | 67%
mdl2 <- lm(Fertility ~ Agriculture+Education+Catholic+Infant.Mortality,swiss)
|=========================================== | 75%
vif(mdl2) Agriculture Education Catholic Infant.Mortality 2.147153 1.816361 1.299916 1.107528
|============================================= | 79%
…
|=============================================== | 83%
1: They are the same. 2: There is no relationship. 3: VIF is the square of standard error inflation.
Selection: 3
|================================================== | 88%
1: We should never exclude anything. 2: We should always exclude it. 3: Excluding it might bias coefficient estimates of regressors with which it is correlated.
Selection: 3
|==================================================== | 92%
1: We should always use uncorrelated regressors. 2: Using converted regressors may make interpretation difficult. 3: Factor analysis takes too much computation.
Selection: 3
|==================================================== | 92%
1: We should always use uncorrelated regressors. 2: Using converted regressors may make interpretation difficult. 3: Factor analysis takes too much computation.
Selection: 2
ANOVA
Finalemte existe una herramienta que nos permite ir viendo progresivamente como evoluciona la inclusión de variables se llama ANOVA, o nested testing models. Se trata de ir añadiendo progresivamente cada uno de los regresores estudiando su impacto.
qownnotes-media-igtMZz
Overfitting and Underfitting (ANOVA)
En el apartado de VIP se demostró cómo incluyendo nuevas variables se incrementa el error estándar de los coeficientes estimados de otras regresores correlacionados. Por tanto, no queremos por no pensar mucho añadir más y más variables dentro del modelo. Por otra parte, omitir varibles puede dar resultado sesgar los coeficientes de los reregresores que están relacionados con los omitidos. A continuación se ven ejemplos del impacto de variables omitidas con el uso de ANOVA para construir representaciones de los datos interpretables y escueto.
…
|========== | 7%
simbias <- function(seed=8765){
# The default seed guarantees a nice histogram. This is the only
# reason that accepting the default, x1c <- simbias(), is required in the lesson.
# The effect will be evident with other seeds as well.
set.seed(seed)
temp <- rnorm(100)
# Point A
x1 <- (temp + rnorm(100))/sqrt(2)
x2 <- (temp + rnorm(100))/sqrt(2)
x3 <- rnorm(100)
# Function to simulate regression of y on 2 variables.
f <- function(k){
# Point B
y <- x1 + x2 + x3 + .3*rnorm(100)
# Point C
c(lm(y ~ x1 + x2)$coef[2],
lm(y ~ x1 + x3)$coef[2])
}
# Point D
sapply(1:150, f)
}
|============== | 11%
1: x2 and x3 2: x1 and x3 3: x1 and x2
Selection: 3
|=================== | 15%
1: 1 2: 0.3 3: 1/sqrt(2)
Selection: 1
|======================== | 19%
x1c <- simbias()
|============================= | 22%
|============================= | 22%
apply(x1c,1,mean) x1 x1 1.034403 1.476944
|================================== | 26%
qownnotes-media-IdkMyQ
# Plot histograms illustrating bias in estimates of a regressor
# coefficient 1) when an uncorrelated regressor is missing and
# 2) when a correlated regressor is missing.
x1hist <- function(x1c){
p1 <- hist(x1c[1,], plot=FALSE)
p2 <- hist(x1c[2,], plot=FALSE)
yrange <- c(0, max(p1$counts, p2$counts))
plot(p1, col=rgb(0,0,1,1/4), xlim=range(x1c), ylim=yrange, xlab="Estimated coefficient of x1",
main="Bias Effect of Omitted Regressor")
plot(p2, col=rgb(1,0,0,1/4), xlim=range(x1c), ylim=yrange, add=TRUE)
legend(1.1, 40, c("Uncorrelated regressor, x3, omitted", "Correlated regressor, x2, omitted"),
fill=c(rgb(0,0,1,1/4), rgb(1,0,0,1/4)))
}
qownnotes-media-OAXbzg
# Illustrate the effect of bogus regressors on residual squared error.
bogus <- function(){
temp <- swiss
# Add 41 columns of random regressors to a copy of the swiss data.
for(n in 1:41){temp[,paste0("random",n)] <- rnorm(nrow(temp))}
# Define a function to compute the deviance of Fertility regressed
# on all regressors up to column n. The function, deviance(model), computes
# the residual sum of squares of the model given as its argument.
f <- function(n){deviance(lm(Fertility ~ ., temp[,1:n]))}
# Apply f to data from n=6, i.e., the legitimate regressors,
# through n=47, i.e., a full complement of bogus regressors.
rss <- sapply(6:47, f)
# Display result.
plot(0:41, rss, xlab="Number of bogus regressors.", ylab="Residual squared error.",
main="Residual Squared Error for Swiss Data\nUsing Irrelevant (Bogus) Regressors",
pch=21, bg='red')
}
…
|================================================ | 37%
fit1 <- lm(Fertility ~ Agriculture,swiss)
|===================================================== | 41%
fit3 <- lm(Fertility ~ Agriculture+Examination+Education,swiss)
|========================================================== | 44%
anova(fit1,fit3)
qownnotes-media-pkQDWT
1: 3102.2 2: 20.968 3: 45
Selection: 2
|=================================================================== | 52%
1: 45 and 43 2: 2 and 3102.2 3: 6283.1 and 3180.9
Selection: 3
|======================================================================== | 56%
deviance(fit3) [1] 3180.925
|============================================================================= | 59%
d <- deviance(fit3)/43
|================================================================================== | 63%
n <- deviance(fit1)-deviance(fit3)
n <- (deviance(fit1) - deviance(fit3))/2
|======================================================================================= | 67%
n/d [1] 20.96783
|=========================================================================================== | 70%
pf(n/d, 2, 43, lower.tail=FALSE) [1] 4.406913e-07
|================================================================================================ | 74%
shapiro.test(fit3$residuals)
Shapiro-Wilk normality test
data: fit3$residuals W = 0.97276, p-value = 0.336
|===================================================================================================== | 78%
qownnotes-media-KXGTlK
|========================================================================================================== | 81%
…
|=============================================================================================================== | 85%
1: Uncorrelated regressors 2: Correlated regressors
Selection: 2
|==================================================================================================================== | 89%
1: False 2: It depends on circumstances. 3: True
Selection: 1
1: True 2: False 3: It depends on circumstances.
Selection: 3
1: False 2: True 3: It depends on circumstances.
Selection: 2
|======================================================================================================================== | 93%
1: Regressors should be tested for normality. 2: Model residuals should be tested for normality. 3: The residuals should be tested for having zero means.
Selection: 3
1: The residuals should be tested for having zero means. 2: Model residuals should be tested for normality. 3: Regressors should be tested for normality.
Selection: 2
TEST3
qownnotes-media-FyoRQF
qownnotes-media-PDMYZn
qownnotes-media-yEfYeY
qownnotes-media-syOIex
qownnotes-media-WugCgL
qownnotes-media-IxTzuQ
GLM (Generalized Linear Models)
Los Modelos Lineales Generalizados es la familia de modelos que incluyen a los modelos lineales, pero incluyen modelos adicionales para los casos en los que son limitados, por ejemplo en respuestas discretas (cómo calculamos el error acumulado) o estrictamente positivas (podríamos usar una transformación logarítmica, o una del arcoseno para los resultados binarios o el logaritmo neperiano para valores sólo positivos), pero con las trasnformaciones muchas veces alteramos los coeficientes de nuestro modelo y son difíciles de interpretar, además de que ya no tenemos los datos en las mismas unidades), por lo que también incluyen regresiones binomiales y binarias así como regresiones de poisson.
Los GLM surgen en 1972 para mejorar las técnicas de transformación. Están formados por 3 componentes:
- Un componente de aleatorieidad proveniente de un modelo de familia exponencial (normal,binomial y poisson) para la respuesta (por ejemplo en los modelos lineales este componente aleatorio eran los errores)
- un componente sistemático via predictor lineal que es la parte que estamos modelando (en modelos lineales son las covariables y los coeficientes del componente lineal)
- Una link functon que conecte la media de la respuesta de la familia exponencial con el preditor lineal. Es fundamental entender aquí que no estaremos transformando la respuesta (las Ys) sino la media de la distribución.
En este curso nos centraremos en las regresiones binomiales o logarítmicas y de poisson, y nos focalizaremos únicamente en las link function más populares.
Ejemplo de regresión lineal * Se asume que \(Y_i \sim N(\mu_i, \sigma^2)\) (la distribución de Gauss es de la familia exponencial.) * Se define el predictor lineal como \(\eta_i = \sum_{k=1}^p X_{ik} \beta_k\). * La link function se define como \(g\) de manera que \(g(\mu) = \eta\). * Para modelos lineales \(g(\mu) = \mu\) de manera que \(\mu_i = \eta_i\) * Esto resulta tener el mismo modelo de verosimilitud likelihood que nuestro error aditivo del modelo linear Gaussiano \[ Y_i = \sum_{k=1}^p X_{ik} \beta_k + \epsilon_{i} \] donde \(\epsilon_i \stackrel{iid}{\sim} N(0, \sigma^2)\)
Regresión de Poisson * Se asume que \(Y_i \sim Poisson(\mu_i)\) de manera que \(E[Y_i] = \mu_i\) donde \(0\leq \mu_i\) * Predictor lineal \(\eta_i = \sum_{k=1}^p X_{ik} \beta_k\) * Link function \(g(\mu) = \eta = \log(\mu)\) * Recordad que \(e^x\) es la inversa de \(\log(x)\) de manera que
\[
\mu_i = e^{\eta_i}
\] Por tanto la verosimilitud es \[
\prod_{i=1}^n (y_i !)^{-1} \mu_i^{y_i}e^{-\mu_i}
\propto \exp\left(\sum_{i=1}^n y_i \eta_i - \sum_{i=1}^n \mu_i\right)
\]
funcionLogistica
- En cada caso, la única manera en la que la verosimilitud depende de los datos es: \[\sum_{i=1}^n y_i \eta_i = \sum_{i=1}^n y_i\sum_{k=1}^p X_{ik} \beta_k = \sum_{k=1}^p \beta_k\sum_{i=1}^n X_{ik} y_i \] Por tanto no necesitamos todos los datos,sólo \(\sum_{i=1}^n X_{ik} y_i\). Esta simplificación es una consecuencia de elegir las llamadas ‘canonical’ link functions.
- (Esto tiene que derivarse). Todos los modelos consiguen su máximo en la raiz de las llamadas ecuaciones normales \[ 0=\sum_{i=1}^n \frac{(Y_i - \mu_i)}{Var(Y_i)}W_i \] donde \(W_i\) son las derivadas de la inversa de la link function.
Una anotación sobre las varianzas, ya que en los modelos lineales la varianza la suponemos constante, mientras que en los casos de bernuilli y de poisson la varianza es variable dependiendo de la observación. En el caso de que las varianzas nos se apeguen lo suficiente a la estructura de varianza definido en el GLM, R proporciona una herramienta llamadas quasi_poisson y quasi que permiten flexibilizar un poco estas restriciones.
- Para el modelo lineal \(Var(Y_i) = \sigma^2\) es constante.
- Para Bernuilli \(Var(Y_i) = \mu_i (1 - \mu_i)\)
- Para Poisson \(Var(Y_i) = \mu_i\).
- En determinados casos, es a menudo relevante tener un modelo de varianza más flexible, incluso si no se corresponde exactamente con la verosimilitud real \[ 0=\sum_{i=1}^n \frac{(Y_i - \mu_i)}{\phi \mu_i (1 - \mu_i ) } W_i ~~~\mbox{and}~~~ 0=\sum_{i=1}^n \frac{(Y_i - \mu_i)}{\phi \mu_i} W_i \]
- Estas son las llamadas ecuaciones normales ‘quasi-likelihood’
Resultados binarios (Bernuilli)
Es muy frecuente encontrarse con resultados que pueden tener únicamente dos posibles valores como vivir o morir, ganar o perder, éxito o fracaso. Estos resultados son llamados binarios, de Bernuilli o 0/1. Una colección de resultados intercambiables binarios para los mismos datos covariables se llaman resultados binomiales. (Los resultados son intercambiables si el orden no importa)
- Bernoulli/binary son usados para modelar salidas con sólo dos posibles resultados
- alive vs dead
- win vs loss
- success vs failure
- disease vs healthy
- binomial outcomes = colección de salidas binarias intercambiables (i.e. tirar una moneda repetidamente) para los mismos datos covariables
- En otras palabras, estamos interesados en conocer el número de \(1\)s predichos vs el número de \(0\)s en lugar de un resultado concreto que sea \(1\) o \(0\)
odds son muy útiles para la construcción de modelos de regresión logistica y son fáciles de interpretar - Vamos a imaginar tirar una moneda con una probabilidad de exito \(p\) + si es cara, ganas \(X\) + si sale cruz, pierdes \(Y\) - Qué valor deberían tener \(X\) e \(Y\) para que el juego sea justo? \[E[earnings]= X p - Y (1 - p) = 0 \Rightarrow \frac{Y}{X} = \frac{p}{1 - p}\] + odds se puede interpretar como “Cuánto estás dispuesto a pagar para una probabilidad \(p\) de ganar 1 dolar?” + si \(p > 0.5\), tienes que pagar más si pierdes que lo que consigues si ganas + si \(p < 0.5\), tienes que pagar menos si pierdes que lo que ganas si aciertas
- odds NO son probabilidades
- odds ratio de 1 = no hay diferencia en odds o 50% - 50%
- \(p = 0.5 \Rightarrow odds = \frac{0.5}{1-0.5} = 1\)
- log odds ratio of 0 = no difference in odds
- \(p = 0.5 \Rightarrow odds = \log\left(\frac{0.5}{1-0.5}\right) = \log(1) = 0\)
- odds ratio < 0.5 or > 2 se llama comunmente un efecto moderado
- Riesgo relativo = tasa de probabilidad en lugar de odds, son a menudo fáciles de interpretar pero más difícil de estimar \[\frac{Pr(W_i | S_i = 10)}{Pr(W_i | S_i = 0)}\]
- Nota: El riesgo relativo a menudo tiene problemas de límite ya el rango de \(\log(p)\) es \((-\infty,~0]\) mientras el rango de logit \(\frac{p}{1-p}\) is \((-\infty,\infty)\)
- Para probabilidades pequeñas el riesgo relativo \(\approx\) Odds Ratio pero no son lo mismo!
Vamos a hacer un ejemplo con predecir la victoria de los Ravens dependiendo del número de puntos que anote. Para eso tenemos un dataset con todos sus resultados.
Ejemplo - Simple Linear Regression
- simple linear regression puedes ser usada para modelar victorias vs derrotas para los Ravens \[ W_i = \beta_0 + \beta_1 S_i + \epsilon_i \]
- \(W_i\) = binary outcome, 1 if a Ravens win, 0 if not
- \(S_i\) = number of points Ravens scored
- \(\beta_0\) = probability of a Ravens win if they score 0 points
- \(\beta_1\) = increase in probability of a Ravens win for each additional point
- \(\epsilon_i\) = residual variation, error
- El valor esperado para el modelo está definido por \[E[W_i | S_i, \beta_0, \beta_1] = \beta_0 + \beta_1 S_i\]
- Sin embargo, el modelo no funcionará bien ya que los resultados a predecir no serán 0 vs 1 sino decimal
- El término de error, \(\epsilon_i\), se presume contínuo y distribuido normal, significando que la predicción probablemente sea decimal
- por tanto, NO es un buen modelo de presunción
# perform linear regression
# load the data
load("ravensData.rda")
head(ravensData)## ravenWinNum ravenWin ravenScore opponentScore
## 1 1 W 24 9
## 2 1 W 38 35
## 3 1 W 28 13
## 4 1 W 34 31
## 5 1 W 44 13
## 6 0 L 23 24
summary(lm(ravenWinNum ~ ravenScore, data = ravensData))##
## Call:
## lm(formula = ravenWinNum ~ ravenScore, data = ravensData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7302 -0.5076 0.1824 0.3215 0.5719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.285032 0.256643 1.111 0.2814
## ravenScore 0.015899 0.009059 1.755 0.0963 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4464 on 18 degrees of freedom
## Multiple R-squared: 0.1461, Adjusted R-squared: 0.09868
## F-statistic: 3.08 on 1 and 18 DF, p-value: 0.09625
- como era de esperar, el modelo produce un pobre ajuste para los datos (\(R^2_{adj} =\) 0.0987)
- Añadiendo un umbral a la salida generada (i.e. si \(\hat W_i < 0.5, \hat W_i = 0\)) y usando el modelo para predecir el resultado sería viable
- sinembargo, los coeficientes para el modelo no son muy interpretables
Recurrimos a la regresión logarítmica
- Probabilidad de victoria de los Ravens quedaría definido por \[Pr(W_i | S_i, \beta_0, \beta_1)\]
- odds queda definido como s \[\frac{Pr(W_i | S_i, \beta_0, \beta_1 )}{1-Pr(W_i | S_i, \beta_0, \beta_1)}\] con rangos de 0 a \(\infty\)
- log odds o logit se define como \[\log\left(\frac{Pr(W_i | S_i, \beta_0, \beta_1 )}{1-Pr(W_i | S_i, \beta_0, \beta_1)}\right)\] yendo de \(-\infty\) a \(\infty\)
- Podemos usar la link function y predictores lineales para construir el modelo de regresion logística \[\begin{aligned} g(\mu_i) & = \log \left(\frac{\mu_i}{1 - \mu_i} \right) = \eta_i\\ (substitute~\mu_i = Pr(W_i | S_i, \beta_0, \beta_1))~g(\mu_i) & = \log\left(\frac{Pr(W_i | S_i, \beta_0, \beta_1)}{1-Pr(W_i | S_i, \beta_0, \beta_1)}\right) = \eta_i \\ (substitute~\eta_i=\beta_0 + \beta_1 S_i) \Rightarrow ~g(\mu_i) & = \log\left(\frac{Pr(W_i | S_i, \beta_0, \beta_1)}{1-Pr(W_i | S_i, \beta_0, \beta_1)}\right) = \beta_0 + \beta_1 S_i\\ \end{aligned} \] que también se puede escribir como \[Pr(W_i | S_i, \beta_0, \beta_1 ) = \frac{\exp(\beta_0 + \beta_1 S_i)}{1 + \exp(\beta_0 + \beta_1 S_i)}\]
- for the model \[\log\left(\frac{Pr(W_i | S_i, \beta_0, \beta_1)}{1-Pr(W_i | S_i, \beta_0, \beta_1)}\right) = \beta_0 + \beta_1 S_i\]
- \(\beta_0\) = log odds of a Ravens win if they score zero points
- \(\beta_1\) = log odds ratio of win probability for each point scored (compared to zero points) \[\beta_1 = \log\left(odds(S_i = S_i+1)\right) - \log\left(odds(S_i = S_i)\right) = \log\left(\frac{odds(S_i = S_i+1)}{odds(S_i = S_i)} \right)\]
- \(\exp(\beta_1)\) = odds ratio of win probability for each point scored (compared to zero points) \[\exp(\beta_1) = \frac{odds(S_i = S_i+1)}{odds(S_i = S_i)}\]
- mediante esta función
manupulatepodemos modificar \(\beta_0\) y \(\beta_1\) y ver como afecta esa modificación la curva de regresión logistica de los datos simulados
# set x values for the points to be plotted
x <- seq(-10, 10, length = 1000)
# "library(manipulate)" is needed to use the manipulate function
manipulate(
# plot the logistic regression curve
plot(x, exp(beta0 + beta1 * x) / (1 + exp(beta0 + beta1 * x)),
type = "l", lwd = 3, frame = FALSE),
# slider for beta1
beta1 = slider(-2, 2, step = .1, initial = 2),
# slider for beta0
beta0 = slider(-2, 2, step = .1, initial = 0)
)Para hacer la regresión logística
# run logistic regression on data
logRegRavens <- glm(ravenWinNum ~ ravenScore, data = ravensData,family="binomial")
# print summary
summary(logRegRavens)##
## Call:
## glm(formula = ravenWinNum ~ ravenScore, family = "binomial",
## data = ravensData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7575 -1.0999 0.5305 0.8060 1.4947
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.68001 1.55412 -1.081 0.28
## ravenScore 0.10658 0.06674 1.597 0.11
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 24.435 on 19 degrees of freedom
## Residual deviance: 20.895 on 18 degrees of freedom
## AIC: 24.895
##
## Number of Fisher Scoring iterations: 5
- Como se puede ver, los coeficientes \(\beta_0\) y \(\beta_1\) son -1.68, 0.107, que son interpretados para ser los log odds ratios
- Podemos convertir los log ratios así como los intervalos de confianza de log a ratios e intervalos de confianza (en las mismas unidades que los datos)
# take e^coefs to find the log ratios
exp(logRegRavens$coeff)## (Intercept) ravenScore
## 0.1863724 1.1124694
# take e^log confidence interval to find the confidence intervals
exp(confint(logRegRavens))## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 0.005674966 3.106384
## ravenScore 0.996229662 1.303304
- Nota: \(\exp(x) \approx 1 + x\) para valores pequeños (cerca de 0) de x, esto puede ser una manera rápida de estimar los coeficientes
- Podemos interpretar la pendiente , \(\beta_1\) como 11.247 % de incremento de probabilidad de ganar por cada punto marcado
- Podemos interpretar el intercept, \(\beta_0\) como 0.186 es el odds de ganar de los Ravens si marcan 0 puntos
- Nota: similar al intercept de un modelo de regresión simple, ek intercept debe ser interpretado con cuidado ya que es un valor extrapolado del modelo y puede que no tenga significado práctico
- para calcular una probabilidad específica de ganar para un número de puntos \[Pr(W_i | S_i, \hat \beta_0, \hat \beta_1) = \frac{\exp(\hat \beta_0 + \hat \beta_1 S_i)}{1 + \exp(\hat \beta_0 + \hat \beta_1 S_i)}\]
- La curva de regresión logistica que resulta es
# plot the logistic regression
plot(ravensData$ravenScore,logRegRavens$fitted,pch=19,col="blue",xlab="Score",ylab="Prob Ravens Win")- ANOVA puede ser llevado a cabo en una regresión logarítmica de un parámetro, donde se analiza el cambio de varianzas con la adición de parámetros al modelo, o multiple nested logistic regression (similar a linear models)
# perform analysis of variance
anova(logRegRavens,test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: ravenWinNum
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 19 24.435
## ravenScore 1 3.5398 18 20.895 0.05991 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
- ANOVA devuelve información sobre el modelo,la link function, la respuesta, así como el análisis de la varianza por anádir términos
Df= change in degrees of freedom- the value 1 refers to adding the
ravenScoreparameter (slope)
- the value 1 refers to adding the
Deviance= measure of goodness of model fit compare to the previous modelResid. Dev= residual deviance for current modelPr(>Chi)= used to evaluate the significance of the added parameter- in this case, the
Deviancevalue of 3.54 is used to find the corresponding p-value from the Chi Squared distribution, which is 0.06- Nota: La distribución Chi Squared con 1 grado de libertad es simplemente la raiz cuadrada de la distribución normal normal , así que el z statistic de 2 corresponde al 95% para distribución normal indica que la desviación de 4 corresponde a approximadamente 5% en la distribución Chi Squared (que es lo que muestra el resultado)
- in this case, the
Ejercicios
The Baltimore Ravens are a team in the American Football League. In post season (championship) play they win about 2/3 of their games. In other words, they win about twice as often as they lose. If I wanted to bet on them, I would have to offer 2-to-1 odds–if they lost I would pay you $2, but if they won you would pay me only $1. That way, in the long run over many bets, we’d both expect to win as much money as we’d lost.
1: 11 to 9 2: 1.22222 to 1 3: Any of these 4: 55 to 45
Selection: 2
|================== | 19%
…
|====================== | 23%
qownnotes-media-ylCMqU
1: 23 2: 30 3: 40 4: 18
Selection: Enter an item from the menu, or 0 to exit Selection: 1
|============================== | 31%
…
|================================= | 35%
ravenData ravenWinNum ravenWin ravenScore 1 1 W 9 2 0 L 13 3 1 W 13 4 1 W 16
qownnotes-media-gDqand
…
|============================================ | 46%
The “best” b0 and b1 are those which maximize the likelihood of the actual win/loss record. Based on the score of a game, b0 and b1 give us a log odds, which we can convert to a probability, p, of a win. We would like p to be high for the scores of winning games, and low for the scores of losses.
…
|================================================ | 50%
We can use R’s glm() function to find the b0 and b1 which maximize the likelihood of our observations. Referring back to the data frame, we want to predict the binary outcomes, ravenWinNum, from the points scored, ravenScore. This corresponds to the formula, ravenWinNum ~ ravenScore, which is the first argument to glm. The second argument, family, describes the outcomes, which in our case are binomial. The third argument is the data, ravenData. Call glm with these parameters and store the result in a variable named mdl.
mdl <- glm(ravenWinNum ~ ravenScore, binomial, ravenData)
|==================================================== | 54%
qownnotes-media-zBwzji
The model is less credible at scores lower than 9. Of course, there is no data in that region; the | Ravens scored at least 9 points in every game. The model gives them a 33% chance of winning if they | score 9 points, which may be reasonable, but it also gives them a 16% chance of winning even if they | score no points! We can use R’s predict() function to see the model’s estimates for lower scores. The | function will take mdl and a data frame of scores as arguments and will return log odds for the give | scores. Call predict(mdl, data.frame(ravenScore=c(0, 3, 6))) and store the result in a variable called | lodds.
lodds <- predict(mdl, data.frame(ravenScore=c(0, 3, 6)))
|=========================================================== | 62%
exp(lodds)/(1+exp(lodds)) 1 2 3 0.1570943 0.2041977 0.2610505
|=============================================================== | 65%
summary(mdl)
Call: glm(formula = ravenWinNum ~ ravenScore, family = binomial, data = ravenData)
Deviance Residuals: Min 1Q Median 3Q Max
-1.7575 -1.0999 0.5305 0.8060 1.4947
Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.68001 1.55412 -1.081 0.28 ravenScore 0.10658 0.06674 1.597 0.11
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 24.435 on 19 degrees of freedom
Residual deviance: 20.895 on 18 degrees of freedom AIC: 24.895
Number of Fisher Scoring iterations: 5
|================================================================== | 69%
exp(confint(mdl)) Waiting for profiling to be done… 2.5 % 97.5 % (Intercept) 0.005674966 3.106384 ravenScore 0.996229662 1.303304
|========================================================================== | 77%
1: 0.005674966 2: 2.5% 3: 0.996229662
Selection: 1
|============================================================================== | 81%
1: They will increase by 30% 2: They will increase slightly 3: They will decrease slightly
Selection: 2
1: They will decrease slightly 2: They will increase by 30% 3: They will increase slightly
Selection: 1
|================================================================================= | 85%
Linear regression minimizes the squared difference between predicted and actual observations, i.e., | minimizes the variance of the residual. If an additional predictor significantly reduces the residual’s | variance, the predictor is deemed important. Deviance extends this idea to generalized linear | regression, using (negative) log likelihoods in place of variance. (For a detailed explanation, see the | slides and lectures.) To see the analysis of deviance for our model, type anova(mdl).
anova(mdl) Analysis of Deviance Table
Model: binomial, link: logit
Response: ravenWinNum
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev
NULL 19 24.435 ravenScore 1 3.5398 18 20.895
|========================================================================================= | 92%
qchisq(0.95, 1) [1] 3.841459
|============================================================================================ | 96%
Resultados de conteo (Poisson)
Las visitas a una web tienden a ocurrir de una manera independiente, una cada vez y en un cierto ratio medio. La distribución de Poisson describe procesos de este tipo. Un proceso de Poisson se caracteriza por un único parámetro, el ratio de ocurrencia esperado, que habitualmente se llama lambda. En nuestro ejemplo, lambda sería las visitas esperadas por día. Por supuesto, si una web se vuelve popular lambde irá creciendo. En otras palabras, lambda dependerá del tiempo. Usaremos la regresión de Poisson para modelar esa dependencia.
- Poisson distribution es un modelo útil para cuentas y tasas
- rate(tasa) = conteos por unidad de tiempo.
- La regresión lineal con transformación es una alternativa.
- Ejemplos de conteos
- Número de casos de gripo en un área
- Número de pacientes que cruzan la puerta de urgencias
- rate data
- Porcentaje de niños que pasan un test
- Porcentaje de entradas en una web desde un país
- grado de radioactividad
- Poisson model examples
- Modelar tasas de entradas en el tráfico de una web
- Aproximación de probabilidades binomiales con \(p\) pequeña y \(n\) grande
- Analizar tablas de contingencia (Conteos tabulados por variables categóricas)
Ejemplo de visitas a una web
Para este ejemplo, el tiempo es siempre un día, así que \(t = 1\), La media de Poisson se interpreta por clicks al día
# laod data
load("gaData.rda")
# convert the dates to proper formats
gaData$julian <- julian(gaData$date)
# plot visits vs dates
plot(gaData$julian,gaData$visits,pch=19,col="darkgrey",xlab="Julian",ylab="Visits")REGRESIÓN LINEAL
El tráfico de la web puede ser modelado linealmente así \[ NH_i = \beta_0 + \beta_1 JD_i + \epsilon_i \] - \(NH_i\) = number of hits to the website - \(JD_i\) = day of the year (Julian day) - \(\beta_0\) = number of hits on Julian day 0 (1970-01-01) - \(\beta_1\) = increase in number of hits per unit day - \(\epsilon_i\) = variation due to everything we didn’t measure * El resultado obtenido es \[ E[NH_i | JD_i, \beta_0, \beta_1] = \beta_0 + \beta_1 JD_i\]
# plot the visits vs dates
plot(gaData$julian,gaData$visits,pch=19,col="darkgrey",xlab="Julian",ylab="Visits")
# perform linear regression
lm1 <- lm(gaData$visits ~ gaData$julian)
# plot regression line
abline(lm1,col="red",lwd=3)TRANSFORMACIÓN LOGARÍTMICA
- Si estamos interesados en el incremento relativo del tráfico, podemos aplicar el logaritmo natural del resultado, de manera que el modelo lineal queda \[ \log(NH_i) = \beta_0 + \beta_1 JD_i + \epsilon_i\]
- \(\log(NH_i)\) = number of hits to the website
- \(JD_i\) = day of the year (Julian day)
- \(\beta_0\) = log number of hits on Julian day 0 (1970-01-01)
- \(\beta_1\) = increase in log number of hits per unit day
- \(\epsilon_i\) = variation due to everything we didn’t measure
- Cuando tomamos el logartimo natural de la salida y se ajusta al modelo de regresión, los coeficientes exponenciados estiman cantidades basadas en las medias geométricas más que en los valores medidos
- \(e^{E[\log(Y)]}\) = geometric mean of \(Y\)
- geometric means are defined as \[e^{\frac{1}{n}\sum_{i=1}^n \log(y_i)} = (\prod_{i=1}^n y_i)^{1/n}\] which is the estimate for the population geometric mean
- as we collect infinite amount of data, \(\prod_{i=1}^n y_i)^{1/n} \to E[\log(Y)]\)
- \(e^{\beta_0}\) = estimated geometric mean hits on day 0
- \(e^{\beta_1}\) = estimated relative increase or decrease in geometric mean hits per day
- Nota: No podemos tomar el logaritmo natural de cero conteos, así que a menudo es necesario aádidr una constante (i.e. 1) para construir el log model
- Añadir la constante cambia ligeramente la inerpretación del coeficiente
- \(e^{\beta_1}\) is now the relative increase or decrease in geometric mean hits + 1 per day
- \(e^{E[\log(Y)]}\) = geometric mean of \(Y\)
- El resultado esperado es \[E[\log(NH_i | JD_i, \beta_0, \beta_1)] = \beta_0 + \beta_1 JD_i \]
round(exp(coef(lm(I(log(gaData$visits + 1)) ~ gaData$julian))), 5)## (Intercept) gaData$julian
## 0.00000 1.00231
Como se puede ver arriba, el incremento diario de clicks es 0.2%
REGRESIÓN DE POISSON
- El modelo de poisson puede definirse como el log de la media \[\log\left(E[NH_i | JD_i, \beta_0, \beta_1]\right) = \beta_0 + \beta_1 JD_i\] or in other form \[E[NH_i | JD_i, \beta_0, \beta_1] = \exp\left(\beta_0 + \beta_1 JD_i\right)\]
- \(NH_i\) = number of hits to the website
- \(JD_i\) = day of the year (Julian day)
- \(\beta_0\) = expected number of hits on Julian day 0 (1970-01-01)
- \(\beta_1\) = expected increase in number of hits per unit day
- Nota: El modelo de poisson difiere de la transformación logarítmica en que los coeficientes son interpretados naturalmente como un valor experado de salida, mientras que en la transformación se interpreta en la escala logarítmica de la salida
- Podemos trasnformar el modelo de Poisson en \[E[NH_i | JD_i, \beta_0, \beta_1] = \exp\left(\beta_0 + \beta_1 JD_i\right) = \exp\left(\beta_0 \right)\exp\left(\beta_1 JD_i\right)\]
- \(\beta_1 = E[NH_i | JD_i+1, \beta_0, \beta_1] - E[NH_i | JD_i, \beta_0, \beta_1]\)
- \(\beta_1\) puede ser por tanto interpretada con el incremento/decremento relativo de clicks por el incremento de 1 día.
glm(outcome~predictor, family = "poisson")= performs Poisson regression
# plot visits vs dates
plot(gaData$julian,gaData$visits,pch=19,col="darkgrey",xlab="Julian",ylab="Visits")
# construct Poisson regression model
glm1 <- glm(gaData$visits ~ gaData$julian,family="poisson")
# plot linear regression line in red
abline(lm1,col="red",lwd=3)
# plot Poisson regression line in
lines(gaData$julian,glm1$fitted,col="blue",lwd=3) …
|========= | 9%
var(rpois(1000, 50)) [1] 50.21732
|============ | 12%
1: The Gauss-Markov BLUE Theorem 2: The Pythagorean Theorem 3: The Central Limit Theorem
Selection: 3
|=============== | 16%
qownnotes-media-cqLKvL
qownnotes-media-dlDSuV
head(hits) date visits simplystats 1 2011-01-01 0 0 2 2011-01-02 0 0 3 2011-01-03 0 0 4 2011-01-04 0 0 5 2011-01-05 0 0 6 2011-01-06 0 0
|=========================== | 28%
class(hits[,‘date’]) [1] “Date”
|================================= | 34%
as.integer(head(hits[,‘date’])) [1] 14975 14976 14977 14978 14979 14980
|==================================== | 38%
mdl <- glm(visits ~ date,poisson, hits)
|======================================= | 41%
…
|========================================== | 44%
summary(mdl)
Call:
glm(formula = visits ~ date, family = poisson, data = hits)
Deviance Residuals:
Min 1Q Median 3Q Max
-5.0466 -1.5908 -0.3198 0.9128 10.6545
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.275e+01 8.130e-01 -40.28 <2e-16 ***
date 2.293e-03 5.266e-05 43.55 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 5150.0 on 730 degrees of freedom
Residual deviance: 3121.6 on 729 degrees of freedom
AIC: 6069.6
Number of Fisher Scoring iterations: 5
|============================================= | 47%
…
|================================================ | 50%
exp(confint(mdl, ‘date’)) Waiting for profiling to be done… 2.5 % 97.5 % 1.002192 1.002399
|=================================================== | 53%
…
|====================================================== | 56%
qownnotes-media-RxDUZg
…
|============================================================ | 62%
which.max(hits[,‘visits’]) [1] 704
|=============================================================== | 66%
hits[704,] date visits simplystats 704 2012-12-04 94 64
|================================================================== | 69%
lambda <- mdl$fitted.values[704] lambda =24.5 | You are amazing!
|===================================================================== | 72%
qpois(.95, lambda) [1] 33
|======================================================================== | 75%
…
|=========================================================================== | 78%
…
|============================================================================== | 81%
qownnotes-media-kIgbwG
mdl2 <- glm(formula = simplystats ~ date, family = poisson, data = hits, offset = log(visits+1))
|================================================================================= | 84%
qpois(.95,mdl2$fitted.values[704]) [1] 47
|==================================================================================== | 88%
1: the square root of lambda. 2: lambda 3: lambda squared
Selection: 3
1: lambda squared 2: the square root of lambda. 3: lambda
Selection: 3
|======================================================================================= | 91%
1: The mean 2: The log of the mean 3: The counts
Selection: 3
1: The log of the mean 2: The counts 3: The mean
Selection: 1
|========================================================================================== | 94%
1: data 2: b0 3: offset 4: family 5: formula
Selection: 3
Ajustes en funciones
Lo que tratamos de hacer ahora es encajar múltiples modelos lineales cuando la curva hace zigzags.
\(Y_i = f(X_i) + \epsilon_i\)
- Se considera el modelo l \[Y_i = \beta_0 + \beta_1 X_i + \sum_{k=1}^d (x_i - \xi_k)_+ \gamma_k + \epsilon_{i}\] donde \((a)_+ = a\) si \(a > 0\) y \(0\) sino y \(\xi_1 \leq ... \leq \xi_d\) con concidos como knot points
- La media \[E[Y_i] = \beta_0 + \beta_1 X_i + \sum_{k=1}^d (x_i - \xi_k)_+ \gamma_k\] es contínua en los knot points
- para \(\xi_k = 5\), el valor esperado oara \(Y_i\) cuando \(x_i\) se aproxima a 5 por la izquierda es \[\begin{aligned} E[Y_i]_{\xi = 5 | left} & = \beta_0 + \beta_1 x_i + (x_i - 5)_+ \gamma_k \\ (since~x_i<5)~ E[Y_i]_{\xi = 5 | left} & = \beta_0 + \beta_1 x_i \\ \end{aligned}\]
- el valor esperado de \(Y_i\) cuando \(x_i\) se aproxima 5 por la derecha es \[\begin{aligned} E[Y_i]_{\xi = 5 | right} & = \beta_0 + \beta_1 x_i + (x_i - 5)_+ \gamma_k \\ (since~x_i>5)~ E[Y_i]_{\xi = 5 | right} & = \beta_0 + \beta_1 x_i + (x_i - 5) \gamma_k\\ (simplify)~ E[Y_i]_{\xi = 5 | right} & = \beta_0 - 5 \gamma_k + (\beta_1 + \gamma_k) x_i \\ \end{aligned}\]
- como se puede ver, el lado derecho es sólo otra linea con diferente intercept (\(\beta_0 - 5 \gamma_k\)) y pendiente (\(\beta_1 + \gamma_k\))
- Mientras \(x\) se aproxima a 5, ambos lados convergen
Example - Fitting Piecewise Linear Function
# simulate data
n <- 500; x <- seq(0, 4 * pi, length = n); y <- sin(x) + rnorm(n, sd = .3)
# define 20 knot points
knots <- seq(0, 8 * pi, length = 20);
# define the ()+ function to only take the values that are positive after the knot pt
splineTerms <- sapply(knots, function(knot) (x > knot) * (x - knot))
# define the predictors as X and spline term
xMat <- cbind(x, splineTerms)
# fit linear models for y vs predictors
yhat <- predict(lm(y ~ xMat))
# plot data points (x, y)
plot(x, y, frame = FALSE, pch = 21, bg = "lightblue")
# plot fitted values
lines(x, yhat, col = "red", lwd = 2)El problema de esta aproximación es que la recta no es contínua ni integrabel en los knot points, si queremos que lo sea deberemos añadir el cuadrado a los términos, y el modelo se transforma en \[Y_i = \beta_0 + \beta_1 X_i + \beta_2 X_i^2 + \sum_{k=1}^d (x_i - \xi_k)_+^2 \gamma_k + \epsilon_{i}\] donde \((a)^2_+ = a^2\) if \(a > 0\) y \(0\)
Añadir términos cúbivos lo haría continuamente integrable en los knot points doblemente,etc..
# define the knot terms in the model
splineTerms <- sapply(knots, function(knot) (x > knot) * (x - knot)^2)
# define the predictors as x, x^2 and knot terms
xMat <- cbind(x, x^2, splineTerms)
# fit linear models for y vs predictors
yhat <- predict(lm(y ~ xMat))
# plot data points (x, y)
plot(x, y, frame = FALSE, pch = 21, bg = "lightblue")
# plot fitted values
lines(x, yhat, col = "red", lwd = 2)TEST4
qownnotes-media-lALHIT
qownnotes-media-TwmpQG
qownnotes-media-CAkQqH
qownnotes-media-EYMoTw
qownnotes-media-keCAhs
qownnotes-media-MiXJKa